Stochastic treatment of finite- A^ effects in mean-field systems and its application to 

the lifetimes of coherent structures 
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A stochastic treatment yielding to the derivation of a general Fokker-Planck equation is presented 
to model the slow convergence towards equilibrium of mean-field systems due to finite- A'^ effects. The 
thermalization process involves notably the disintegration of coherent structures that may sustain 
out-of-equilibrium quasistationary states. The time evolution of the fraction of particles remaining 
close to a mean-field potential trough is analytically computed. This indicator enables to estimate 
the lifetime of coherent structures and thermalization timescale in mean-field systems. 
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Many physical systems may be considered as isolated 
assemblies of N bodies interacting via long-range pair 
interactions. This is the case for systems ranging from 
charged particles interacting via Coulomb interaction to 
self-gravitating massive objects like globular clusters or 
stars in galaxies and this may even include suitably pre- 
pared Bose-Einstein condensates 'l| in a close future. 
The physically relevant issue of the dynamics of those sys- 
tems in the large- A^ limit forms the subject of kinetic the- 
ory. Long-range systems are prone to collective behavior 
that may be largely dominating before binary collisional 
effects set on. This hierarchy between collective and col- 
lisional behavior is responsible for the unusual properties 
of the relaxation process towards equilibrium as well as 
for the richness and complexity of the physics of long- 
range systems. These are motivations for the present con- 
siderable interest raised by long-range systems in various 
fields such as plasma physics, astrophysics [2], statistical 
physics [3J or applied mathematics. 

Collective behavior of long-range systems as well as the 
intricacies of the relationships between their dynamics, 
kinetic theory and equilibrium statistical properties may 
be more conveniently unveiled through models that are 
already of mean-field type for finite N. These are Hamil- 
tonian models describing e.g. wave-particle interaction 
[J] , which is an ubiquitous phenomenon in hot and dilute 
plasmas, or the all-to-all coupling of A^ bodies in long- 
range interactions of the kind of gravitation in a com- 
pact space [si |6| . Despite their relative simplicity, such 
models develop a rich long-range phenomenology. This 
includes in particular the emergence of quasistationary 
states (QSSs) having lifetimes diverging with A^, during 
which the time average of macroscopic quantities, such 
as the temperature or the modulus of mean-fields, differs 
from their equilibrium statistical mechanics ensemble av- 
erages. These QSSs may be connected to the existence 
of coherent structures that may be viewed as long-lived 
phase space patterns related to locally insufficient mix- 
ing properties [7[. Consequently, the relaxation to equi- 
librium should accompany the disintegration of coherent 



structures. 

The Letter is organized as follows: first, a stochastic 
treatment of finite A^-effects in mean field systems will 
be proposed leading to the establishment of a Fokker- 
Planck equation. In order to test this model, we then 
shall consider an Hamiltonian model of A^ particles in 
self-consistent interaction via a cosine potential. Start- 
ing from configurations where 0{N) particles are trapped 
into their self-potential well, an analytic expression giv- 
ing the fraction of the particles that remain trapped as 
a function of time will be successfully tested against nu- 
merical results. The relevance of this indicator to the 
thermalization issue will be shortly discussed. 

Consider A^ particles evolving in the phase space S^ x 
M.^ under the dynamics deriving from the Hamiltonian 



N 



^-E 



i=l 



Pi 

2 



1 

2N 



N 



E E ^" ^°^ [^"(^J ^ ^*)] ' (^) 



n— 1 i,j — 1 



where qi £ Sl is the position of particle i on the cir- 
cle Sl = R/i, Pi its conjugate momentum, and where 
only the first s long-range components with wave num- 
bers kn = 27Tn/L, for 1 < n < s, are retained in the 
potential term. When y„ ex n~^, model (H]) amounts to 
the long-range truncation of the one-dimensional Newto- 
nian potential with space-periodic boundary conditions, 
that describes Coulomb or gravitational interaction de- 
pending on the potential sign. Various systems covered 
by H]) have been discussed in Ref . [8| . An extension to a 
spatial dimension d > 1 should not be a conceptual prob- 
lem. Introducing the set of collective observables {M„} 
through 
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M„ = -r^ E (cos(fc„(jj),sin(fc„(jj)) = M„(cos(?!)„,sin0„) 

(2) 
yields the equation of motion of any particle i as 

s 

Qi + ^ knVnMn sin(fc„gj - (/>„) = 0. (3) 

n=l 



Therefore, the collective variables M„ behave as mean 
fields that, as well as the phases </)„, depend on time 
through the self-consistency relations © . 

For smooth potentials like in ([T|), the convergence of 
the finite- A'^ dynamics to Vlasov equation is rigorously 
proved on arbitrary finite-time intervals 9] . Vlasov equa- 
tion being time reversible, its solution f(j),q,t) cannot 
approach an equilibrium [101] . yet macroscopic quantities, 
involving phase space integrals of /, such as the mean- 
fields, can converge to stationary values associated to 
QSSs. In the realistic finite-A^ Hamiltonian framework, 
finite- A'^ effects will in the long term induce the thermal- 
ization process and disintegration of coherent structures 
possibly sustaining those QSSs. 

Modelling this process, we assume that the system ([T]) 
is trapped in a QSS, such that one can write the mean 
fields as M„{t) = M^(t) + SMnit) with SMn < A/O, 
where SMn varies on a time scale which is very much 
smaller than the characteristic time scale of M°, the lat- 
ter being comparable to a local average. We now make 
the hypothesis that during the QSS regime, the fluctu- 
ations around the local mean value have a variance de- 
creasing as N^^. The central idea behind this is to re- 
place the deterministic but yet very chaotic fluctuations 
of SMn by stochastic processes, whose variances are suit- 
ably chosen. Numerical observations support this model- 
ing. For instance, as shown in Fig.[T]for some special case, 
the mean fields clearly exhibit two different timescales, 
in agreement with the decomposition suggested earlier. 
Moreover, both histograms and the autocorrelation func- 
tion shown in Fig. [T] suggest that one can model the fluc- 
tuations SMn by a Gaussian (white) noise. 



the stochastic version of the equations of motion ([3]) reads 
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FIG. 1: (Left) Comparison between the numerically com- 
puted fluctuations of Mi, M2 and M3 around their local av- 
erage (Mn) with a Gaussian fit, for system 11} with s = 3 
and V„ = n starting from an arbitrary initial condition. The 
inset shows the evolution of Mi, M2 and M3 with respect 
to time. (Right) Autocorrelation of Mi after subtracting its 
local average. 



Replacing SMn by the noise ^„ such that {£,n)—0, and 
{Cn{t)^n{t')) cx N-^S(t-t'), the Fokker-Planck equation 
(FPE) associated to the Langevin equations coming from 
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This equation may be interpreted as a Vlasov equation 
supplemented with a r.h.s. of order N~^, consistently 
with the argument presented in Ref. |3| , coming here not 
from binary collisions but from the fluctuations of the 
mean flelds. This differs from other FPEs derived in 
mean-field systems in other places [ll|, \l^. Moreover, 
this equation was derived without any need to invoke 
dissipation (see e.g. [131]). 

In what follows, we shall consider the case of a sin- 
gle resonance (s = 1) in which coherent structures may 
survive for long times close to the potential trough (see 
e.g. Fig. [2]). The FPE ^ may be further simplified by 
looking for solutions in separate variables q and p writing 
f{p,q,t) = g{q,t)f(p,t). Assuming that g is even in the 
wave frame and that J„ g[q, t)dq is constant, one obtains 
a simple diffusion equation 
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'{kq~(j))g{q,t)dq/ ^Q g{q,t)dq. 
Numerical evidence [3| supports the fact that p and q may 
be treated as separate variables and that q may be con- 
sidered as a fast variable compared to p, meaning that the 
distribution function in q approaches much more quickly 
its Boltzmann-Gibbs shape than the p one consistently 
with basic dimensional arguments 13[ . The forthcoming 
numerical tests will show that the average of sin^ {kq~<p) 
may be effectively replaced by its ensemble average. 

From now on, in order to simplify expressions, we shall 
put V — 1 and L — 27r, which gives k = 1. This amounts 
to the well-known Hamiltonian Mean Field model [3|. 
Figure [5] shows numerical results obtained starting from 
a monokinetic beam [6| . The upper panel shows the one- 
particle phase space plots at three different stages of the 
evolution. Deep inside the mean-field potential trough, 
particles move almost regularly forming a clear coherent 
structure pattern that progressively dissipates. It is in- 
teresting to note the similarity of these figures with the 
phase space plots for the 1-D finite- A^ cold dark matter 
simulations of Ref. [IJ] . The lower panel shows the evolu- 
tion of the mean-field for four different numbers of parti- 
cles, ranging from A^ = 10^ to A^ = 10^. The dashed line 
marks its equilibrium ensemble average. It is clear from 







FIG. 2: (Upper panel) Phase-space snapshots at increasing 
times (from left to right, i = 6, i = 100 and t = 2000) 
for an initial cold-beam condition of energy U = 0.5. The 
simulation was performed using A'' — 10* particles. The in- 
stantaneous separatrix is plotted in bold red. (Lower panel) 
Time evolution of the magnetization M and fraction ni/N 
for N = 10^ (green), 2.10^ (red), 5.10^ black and 10'' (blue) 
particles. 



the figure that the convergence towards equilibrium slows 
down as N increases. Also displayed is the evolution of 
ni/N , that is the fraction of particles initially trapped 
within the separatrices and which remain inside them up 
to time t. During a short transient regime (not shown), 
about twenty percent of the particles escape the mean- 
field resonance. It is only after this lapse of time that 
the regime becomes diffusive. Most importantly. Fig. [2] 
gives an evidence that the cancelation time of ni/N may 
be used as a marker of thermalization as it qualitatively 
coincides with the time where the mean-field begins to 
fiuctuate around its equilibrium value in a stable way. 
This is not surprising since, at the time when ni/N van- 
ishes, the coherent structure has been completely disin- 
tegrated. The dependence of M on A^ becomes much 
clearer when increasing iV to values of the order of 10^, 
but the numerical cost was too important for us to per- 
form a full simulation leading to a vanishing rii/N . How- 
ever, the early-time behaviour of ni/N was still correctly 
described by our model. 

The equations of the separatrices are given by ps {t) — 
±2^jM{t)cos{q/2). Let us put A = 2\/M" and consider 
it to be almost constant. Figure [5] is a motivation to 
answer the following question: What is the number of 
particles ni{t), having initially momenta comprised be- 
tween —A and A, that remain in this domain up to time 
t ? Among such particles are the particles forming the 
core of coherent structures that slow down the mixing 
process. However, because of the parametric resonance 



induced by the mean fields fluctuations, these particles 
will eventually escape, and we shall now estimate the 
characteristic time needed by the system to evacuate a 
fraction 1 — S of the A^o particles initially contained in 
the band of momenta [—A; A]. Using the linearity of the 
diffusion equation ([5]) , one can focus on the contribution 
of particles whose momenta remain up to time t within 
the band [—A; A] and solve ([5]) by imposing the cancela- 
tion of / at p = ±A. Finally, integrating over p yields 
the solution 
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with T) given by Eq. ^ and 
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As confirmed by the numerical simulations shown in 
Fig. [31 the dynamics arising from Eq. ([5]) proves to cor- 
rectly depict the escape process through the validation 
of Eqs. ([2l)-([8]). In these numerical simulations, par- 
ticles were initially distributed according to a so-called 
waterbag distribution 

fo{p,q)^e{Ap-\p\)e{Aq-\q\)/{AApAq) (9) 

where Q stands for the Heaviside step function. It is 
interesting to note that when starting from these initial 
waterbag conditions, the phase-space distribution even- 
tually exhibits a core-halo structure, which has been 
recently investigated in [18[. Once again, it is visible 





FIG. 3: (Left) Snapshot of the one-particle phase space at 
i = 12 for A'' — 5000 particles initially distributed in a wa- 
terbag configuration ([9]) with Ap = 0.848 and Aq = 2.16. 
(Right) Time evolution of the mean-field M (red curve) and 
of Tii/N (black curve). The blue dashed curve is the analytic 
expression for ni/N as deduced from Eqs. ©-((S]!. The initial 
time has been chosen at time i = 10 and A'^o = 0.95A''. 

on Fig. [3] that, at the time when all the particles that 
where initially in the momentum band [—A; A] have at 
least once escaped this domain, the system has seemingly 
reached its thermal equilibrium. In order to test the es- 
cape model given by Eqs. (©-([S]), one needs to know 
the diffusion coefficient V given by Eq. ([B]). As already 



discussed, sin^(g) may be estimated from (sin^ g) = 
h W {M)J /[/3 (M)^ /o i/3 (M) J]. A priori (f^) has to be 
determined froni numerical simulations since the system 
is not at equilibrium. However, in the cases that were 
considered, the numerically computed variance (^^) was 
almost indistinguishable from its canonical value given 
by 
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13 (M), sat- 
0. Numeri- 



where V'(^) = ^^7(2/3) — log/o(w) and v" 

isfies the self-consistency equation dyiplv 

cally, this gives (^^) « 0.43/iV, which is consistent with 

the fit obtained from numerical simulations. As shown in 

Fig- El the agreement between the numerically computed 

time evolution of ng/N and its analytic modeling ([II)-® 

is quite satisfactory. 

Let us finally estimate the time needed to destroy the 
inner coherent structure by the means of Eq. ([7]). Con- 
sidering that No = 0{N) and V = O(iV-i), a rough es- 
timate of the time ts{N) needed for ng/N to reach down 
a sufficiently small fraction 5 gives 



n{N) (X -NlogS. 



(11) 



When V depends on time, Eq. ([Tl]) follows from the mean 
value theorem. This first result recovers the linear N 
scaling found in the abundant literature of long-range 
interacting systems [5|, ll5[ ll6| , where the numerical ev- 
idence is extracted from thresholds equivalent to the S 
criterion imposed here. One does not expect the contin- 
uum approach behind Eq. ([7]) to remain valid for vanish- 
ingly small values of 5. However, going up to the limit of 
validity of this model, one may infer that the sweeping 
of phase space is sufficient to reach a complete thermal- 
ization at a time tqss when ng = 0(1). Then Eq. (fTTj) 
would give the maximal scaling tqss oc NlogN. This 



corresponds to the scaling recently suggested in Ref. [17 1 
for the s = 1 case. Eq. pT|) predicts a linear behavior 
with respect to A'^, which we found to be correct over the 
whole range of values of N studied here, independently 
from the chosen threshold 6. We also checked the scaling 
with the latter parameter. Figure |4] shows a very good 
agreement between Eq. (ITT|) and the numerical simula- 
tions. 

These scalings contrast with the numerically obtained 
A^ ^-^ scaling for the QSS lifetime starting from two special 
initial conditions [19|, |20| . These cases are however not 
in contradiction with the results presented here since, 
even this is less obvious for [l9[ , they both correspond to 
QSS about a vanishing mean-field, a case that is excluded 
from the present framework since the phase would be 
no longer defined. In the intermediate cases, where the 
QSS magnetization is clearly above zero, but yet far from 
the equilibrium expectation, this method provides a good 
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FIG. 4: (Left) Plot of ts{N) in log- log scale with respect 
to A'^. As predicted by Eq. (|11[) . the behavior is linear with 
A^. (Right) Plot of the numerically measured rg with respect 
to log 5. The behavior is linear over a wide range of 5. As 
expected, the measured times for very low threshold values 
are lower than the logarithmic prediction of Eq. (|ll|l . since 
the latter is obtained in the continuous limit. 



estimation of the time needed to destroy the coherent 
structures, but fails to predict the QSS lifetime, since 
the effective model does not capture the average growth 
of the separatrix with time. 

The present framework and results are expected to 
be easily transposable to wave-particle models in which 
finite- A^ effects eventually drive the system towards equi- 



librium in contradiction with the Vlasov approach |21 



This discreteness effect may be more than a numerical 
concern for simulations since some physical effects 22 1 
cannot be explained in the Vlasov limit. 
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